% Copyright (C) 2019-2023 Benjamin Born, Francesco D'Ascanio, Gernot J. Mueller, Johannes Pfeifer
%
% This is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% It is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% For a copy of the GNU General Public License,
% see <http://www.gnu.org/licenses/>.

% Run LP test conditioning on slack

function run_state_test_slack(sign_iter)

if sign_iter==1
    sign_dummy = -1;
elseif sign_iter==2
    sign_dummy = 1;
else
    error('Case not implemented')
end

addpath('../Auxiliary_Files')

lag_number=4;
max_irf_horizon = 8;
time_fixed_effects_dummy=1;
country_fixed_effects_dummy=1;

add_residual_dummy=1; %adds residuals from previous step to regression to increase efficiency;


impulse_name = 'PF_shock_G';
special_name = 'pf_fc';

%% Construct sample

% Load data arrays for regression
load('../Figure 3+6 + Appendix C Descriptives/dataset_BDMP.mat')
temp=load('../Sign_restriction_first_stage_shocks_euro_area/results/structural_shocks_point.mat');

data_array_for_regression_stacked_by_variable = cat(3,data_array_for_regression_stacked_by_variable,temp.structural_shocks.BP_G,temp.structural_shocks.BP_T,temp.structural_shocks.PF_G,temp.structural_shocks.PF_T);
pos.BP_shock_G_CK=length(fieldnames(pos))+1;
Header{1,pos.BP_shock_G_CK}='BP_shock G based on Caldara/Kamps';
pos.BP_shock_T_CK=length(fieldnames(pos))+1;
Header{1,pos.BP_shock_T_CK}='BP_shock T based on Caldara/Kamps';
pos.PF_shock_G=length(fieldnames(pos))+1;
Header{1,pos.PF_shock_G}='Sign restriction shock G based on Caldara/Kamps';
pos.PF_shock_T=length(fieldnames(pos))+1;
Header{1,pos.PF_shock_T}='Sign restriction shock T based on Caldara/Kamps';

% check whether positions are unique
pos_numbers=cell2mat(struct2cell(pos));
unique_pos_numbers=unique(pos_numbers);
if ~isequal(pos_numbers,unique_pos_numbers)
    pos_numbers(~ismember(pos_numbers,unique_pos_numbers))
    error('The position numbers are wrong')
end
% Set dependent variables and regressors
regressor_var_names={'g_real_demeaned_linearly_detrended', 'y_real_demeaned_linearly_detrended', 'Effective_FX_Real_Intra_Euro_CPI_log', 'tax_revenue_real_demeaned_detrended'};
dependent_var_names={'g_real_demeaned_linearly_detrended','y_real_demeaned_linearly_detrended','Effective_FX_Real_Intra_Euro_CPI_log', 'tax_revenue_real_demeaned_detrended'};
special_name = ['real_fx_intraeuro_',special_name];
dependent_var_plottitles={'Government consumption','GDP','Real Effective Exchange Rate', 'Tax revenues'};
dependent_var_ylabels={'percent','percent','percent', 'percent'};

if lag_number<1
    error('Lag number must be strictly positive')
end

% Warning: regressor_var_names set to {} below!!

indicator_save_name='ecdf_rec';

%euro countries only
[row,col] = find(data_array_for_regression_stacked_by_variable(:,:,pos.euro_country)==0);
condition=' find(data_array_for_regression_stacked_by_variable(:,:,pos.euro_country)==0);';
split_save_name='eur_slack_interact';
Figure_name_split='Euro Countries';
        
for ii=1:length(row)
    data_array_for_regression_stacked_by_variable(row(ii),col(ii),3:end)=NaN;
end


%% Shock Variable
fe_shocks_temp = data_array_for_regression_stacked_by_variable(:,:,pos.(impulse_name));

fe_shocks_negative = zeros(size(fe_shocks_temp));
fe_shocks_negative(fe_shocks_temp<0 | isnan(fe_shocks_temp)) = fe_shocks_temp(fe_shocks_temp<0 | isnan(fe_shocks_temp));
fe_shocks_positive = zeros(size(fe_shocks_temp));
fe_shocks_positive(fe_shocks_temp>0 | isnan(fe_shocks_temp)) = fe_shocks_temp(fe_shocks_temp>0 | isnan(fe_shocks_temp));
if any(any(fe_shocks_negative>0)) || any(any(fe_shocks_positive<0 ))
    error('Sign is wrong')
end

slack_high = data_array_for_regression_stacked_by_variable(:,:,pos.Indicator_unemployment_rate_above_percentile_50);
slack_low = 1 - slack_high;

if sign_dummy==-1
    data_array_for_regression_impulse_only=cat(3,fe_shocks_negative .* slack_high, fe_shocks_negative .* slack_low,...
        fe_shocks_positive.*slack_high, fe_shocks_positive.*slack_low);
elseif sign_dummy==1
    data_array_for_regression_impulse_only=cat(3,fe_shocks_positive .* slack_high, fe_shocks_positive .* slack_low ,...
        fe_shocks_negative.*slack_high, fe_shocks_negative.*slack_low);
elseif sign_dummy==0
        data_array_for_regression_impulse_only=cat(3,fe_shocks_temp);
else
    error('undefined case for sign_dummy')
end

data_array_for_regression=data_array_for_regression_impulse_only;

%% Lagged regressors
for var_iter=1:length(regressor_var_names)
    data_array_control=data_array_for_regression_stacked_by_variable(:,:,pos.([regressor_var_names{var_iter},'_lag_1']):pos.([regressor_var_names{var_iter},'_lag_',num2str(lag_number)]));
    data_array_for_regression=cat(3,data_array_for_regression, repmat(slack_high,[1,1,4]).*data_array_control, repmat(slack_low,[1,1,4]).*data_array_control);
end

%% Run regression

theta_mat=NaN(1+max_irf_horizon,1,length(dependent_var_names));
se_mat=NaN(1+max_irf_horizon,1,length(dependent_var_names));
coeff_of_interest = size(data_array_for_regression_impulse_only, 3);
theta_large = NaN(1+max_irf_horizon, coeff_of_interest, length(dependent_var_names));
theta_table = NaN(8, 10);

for var_iter=1:length(dependent_var_names)
    for horizon_iter = 0:max_irf_horizon
        %construct matrices without NaN
        if horizon_iter==0
            dependent_variable=squeeze(data_array_for_regression_stacked_by_variable(:,:,pos.(dependent_var_names{var_iter})));
            yhat_full=NaN(size(dependent_variable));
            [dependent_variable_for_run,data_array_for_regression_for_run,time_indices_non_NaN,country_indices_non_NaN,reg_matrices]=...
                create_regression_matrices_no_NaN_indicator(dependent_variable,data_array_for_regression,data_array_for_regression_stacked_by_variable,pos,country_indicator_names_mapping,time_fixed_effects_dummy,country_fixed_effects_dummy,slack_high);
            %run regression
            [theta,stdDK,~,CovDK,yhat] = HszDk5cPs(dependent_variable_for_run,ones(size(dependent_variable_for_run,1),1),data_array_for_regression_for_run,1,7,1);
            
            if add_residual_dummy && ~strcmp(impulse_name,dependent_var_names(var_iter)) %when regression G on itself, residuals are 0
                yhat_full(time_indices_non_NaN,country_indices_non_NaN)=yhat;
                resids=dependent_variable-yhat_full;
            else
                resids=[];
            end
        else
            dependent_variable=squeeze(data_array_for_regression_stacked_by_variable(:,:,pos.([dependent_var_names{var_iter},'_lead_',num2str(horizon_iter)])));
            yhat_full=NaN(size(dependent_variable));
            

            [dependent_variable_for_run,data_array_for_regression_for_run,time_indices_non_NaN,country_indices_non_NaN,reg_matrices]=...
                create_regression_matrices_no_NaN_indicator(dependent_variable,cat(3,data_array_for_regression,resids),data_array_for_regression_stacked_by_variable,pos,country_indicator_names_mapping,time_fixed_effects_dummy,country_fixed_effects_dummy,slack_high);
            %run regression
            [theta,stdDK,~,CovDK,yhat] = HszDk5cPs(dependent_variable_for_run,ones(size(dependent_variable_for_run,1),1),data_array_for_regression_for_run,1,7,1);
            if add_residual_dummy
                yhat_full(time_indices_non_NaN,country_indices_non_NaN)=yhat;
                resids=dependent_variable-yhat_full;
            else
                resids=[];
            end
        end
        
        %save coefficients
        theta_mat(1+horizon_iter,:,var_iter) = theta(1,1);
        se_mat(1+horizon_iter,:,var_iter)  = stdDK(1,1)';
        Rmat = zeros(1,length(theta));
        Rmat(1,1:2) = [1,-1];
        Wstat = (Rmat*theta)'/(Rmat*CovDK*Rmat')*(Rmat*theta);
        test_stat_mat(1+horizon_iter, var_iter) = Wstat;
        test_p_mat(1+horizon_iter, var_iter) = 1 - chi2cdf(Wstat,1);

        if horizon_iter<5 && var_iter==2
            theta_table(1,1+horizon_iter) = -100 * theta(1,1);
            theta_table(2,1+horizon_iter) = 100 * stdDK(1,1);
            theta_table(3,1+horizon_iter) = (1-normcdf(abs(theta(1,1))./stdDK(1,1)))*2;

            theta_table(4,1+horizon_iter) = 100 * theta(2,1);
            theta_table(5,1+horizon_iter) = 100 * stdDK(2,1);
            theta_table(6,1+horizon_iter) = (1-normcdf(abs(theta(2,1))./stdDK(2,1)))*2;

            theta_table(7,1+horizon_iter) = 100 * (theta(1,1) - theta(2,1));
            theta_table(8,1+horizon_iter) = test_p_mat(1+horizon_iter, var_iter);

        elseif horizon_iter<5 && var_iter==3
            theta_table(1,1+horizon_iter+5) = 100 * theta(1,1);
            theta_table(2,1+horizon_iter+5) = 100 * stdDK(1,1);
            theta_table(3,1+horizon_iter+5) = (1-normcdf(abs(theta(1,1))./stdDK(1,1)))*2;

            theta_table(4,1+horizon_iter+5) = -100 * theta(2,1);
            theta_table(5,1+horizon_iter+5) = 100 * stdDK(2,1);
            theta_table(6,1+horizon_iter+5) = (1-normcdf(abs(theta(2,1))./stdDK(2,1)))*2;

                theta_table(7,1+horizon_iter+5) = -100 * (theta(1,1) - theta(2,1)); %minus to account for definition of FX in Eurostat data, which is the opposite of the one in the paper
            theta_table(8,1+horizon_iter+5) = test_p_mat(1+horizon_iter, var_iter);
        end
        if horizon_iter==5 && var_iter==3
            if sign_dummy==1
                title_name='Differences in response coefficients across states and shock signs: row 2';
                headers_string=char('\psi_{h}^{+}|u^h - \psi_{h}^{-}|u^l','0','1','2','3','4');
            else
                title_name='Differences in response coefficients across states and shock signs: row 3';
                headers_string=char('\psi_{h}^{-}|u^h - \psi_{h}^{+}|u^l','0','1','2','3','4');
            end
            labels_string=char('Output','p-val','REER','p-val');
                dyntable(title_name,headers_string,labels_string,[theta_table(7:8,1:5); theta_table(7:8,6:end)],size(labels_string,2)+2,5,2)
        end

        theta_large(1+horizon_iter,:,var_iter) = theta(1:coeff_of_interest,1);
    end
end